Time-dependent Gutzwiller approximation for the Hubbard model 
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We develop a time-dependent Gutzwiller approximation (GA) for the Hubbard model analogous to 
the time-dependent Hartree-Fock (HF) method. The formalism incorporates ground state correla- 
tions of the random phase approximation (RPA) type beyond the GA. Static quantities like ground 
state energy and double occupancy are in excellent agreement with exact results in one dimension 
up to moderate coupling and in two dimensions for all couplings. We find a substantial improve- 
ment over traditional GA and HF+RPA treatments. Dynamical correlation functions can be easily 
computed and are also substantially better than HF+RPA ones and obey well behaved sum rules. 



The Gutzwiller (GW) trial wave function |jj is prob- 
ably the most popular variational approach to the Hub- 
bard model which incorporates correlation effects beyond 
the Hartree-Fock (HF) approximation. Since the Hub- 
bard model describes the competition between hopping 
and correlation induced localization of the charge carri- 
ers, the idea is to apply a projector to a Slater determi- 
nant (SD) which reduces the number of doubly occupied 
sites. The optimum double occupancy probability is de- 
termined variationally. 

Similar to HF best results are obtained if one allows 
for unrestricted charge and spin distributions which are 
determined also variationally. For example for the half- 
filled Hubbard model a SD with long range antiferromag- 
netic order is favoured [Q. 

A formal diagrammatic solution of the GW variational 
problem has been given by Metzner and Vollhardt ||[|, 
however, for the most part of practical purposes one ap- 
proximates the corresponding expectation values using 
the so-called Gutzwiller approximation (GA) 

The GA can be derived using a variety of methods 
In particular it is recovered at the mean- field 
level (saddle-point) of the four-slave boson functional in- 
tegral method introduced by Kotliar and Ruckenstein 
(KR) H . The latter offers the possibility of going be- 
yond the Gutzwiller result as for example the inclusion 
of transversal spin degrees of freedom In addition 
it provides a scheme to include fluctuations beyond the 
mean-field (MF) solution. Expansions around the slave- 
boson saddle point have been performed for homogeneous 
systems in Ref. H|| in order to calculate correlation func- 
tions in the charge and longitudinal spin channels. How- 
ever, the expansion of the KR hopping factor z SB is a 
highly nontrivial task both with respect to the proper 
normal ordering of the bosons and also with respect to 
the correct continuum limit of the functional integral 

The complexity of the expansions around the slave bo- 
son saddle point have severely hampered practical com- 
putations of dynamical quantities within this formalism. 



One of the few successful attempts is the computation 
of the optical conductivity in the paramagnetic state in 
Ref. [ p~2|JT^] . Remarkably although the starting SD de- 
scribes a paramagnetic system, spectral weight on the 
Hubbard bands appears as an effect of fluctuations. As 
far as we know this approach has not been pursued in 
broken symmetry states due to technical difficulties, in- 
cluding the fact that the KR choice for the z SB hopping 
factor does not lead to controlled sum rules | Q . 

In this work we introduce a simple scheme to com- 
pute fluctuation corrections around the GA to dynami- 
cal and static correlation functions and the ground state 
energy. The method can be viewed as a time-dependent 
GA in the same way as the random phase approxima- 
tion (RPA) method on top of a HF solution (HF+RPA) 
can be viewed as time-dependent HF approximation in 
the limit of small amplitude oscillations [ p~5|Jl~6| . For this 
reason we label the method as GA+RPA. It is also a gen- 
eralization of the method of Ref. jlTj in order to describe 
the low temperature Fermi liquid regime. The method 
incorporates ground state correlations beyond the ones 
of the Gutzwiller type just as HF+RPA takes into ac- 
count ground state correlations not present in the HF 
wave function. 

The GA+RPA ground state energy of the one-band 
Hubbard model is in excellent agreement with exact re- 
sults up to moderate coupling in one dimension (Id) and 
for all couplings in a 2d system (Fig. |l|) . The optical con- 
ductivity of a Hubbard chain is in much better agreement 
with numerical results than HF+RPA (Fig. ||). In addi- 
tion sum rules are well behaved in the HF+RPA sense 



We consider the one-band Hubbard model 



H ~ ^""^ tij C i,<j C j,<? 



n ht n hl 



(1) 



where Cj i<T destroys an electron with spin a at site i, and 
n-i,a = cj a Ci. a . U is the on-site Hubbard repulsion and iy 
denotes the transfer parameter between sites i and j. In 
the numerical computations below we take only nearest 
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neighbour matrix elements U 



-t to be non-zero. 



Our starting point is an energy functional E[p,D] of 
the GA type §. Here p is the density matrix of an 
associated Slater determinant \SD >, i.e. Pia,jcr> = 
(SD\cl a Cj a '\SD) and D is a vector of the GA double 
occupancy parameters Di at site i. In order to consider 
arbitrary fluctuations the charge and spin distribution of 
p and the distribution of D should be completely unre- 
stricted. For simplicity we consider only solutions where 
the associated SD is an eigenstate of the z-component of 
the total spin operator (piaja' = S a .a-' Pija-)- 

E[p, D] can be obtained by exploiting the equivalence 
between the KR saddle point solution and the GA Q . It 
is given by: 



E[p, D} = ) hjZivZjoPijo + Uy^Dj 



(2) 



and 



V(l - pa + Di)( Piia - Di) + y/Di( Pii - a - Di) 



\J Piicri\ ~ Piicr) 



(3) 



with pa = pucr- The stationary solution p^°',D^ 
is determined by minimizing the energy functional with 
respect to p and D. 

The variation with respect to the density matrix has to 
be constrained to the subspace of Slater determinants by 
imposing the projector condition p 2 — p Within 
this subspace we now consider small time-dependent am- 
plitude fluctuations of the density matrix p(t). We add 
a weak time-dependent field to Eq. (0) of the form: 

F ( f ) = Y,icr,ja'(fi^"' e ^ iUtc L c j^ +h.c). This produces 
small amplitudes oscillations 8 pit) around the stationary 
density i.e. Spit) = pit) — p 



(0) 



We assume that at each instant of time the double oc- 
cupancy parameter is at the minimum of the energy func- 
tional compatible with the corresponding p(t); i.e. the 
double occupancy parameters {D} adjust antiadiabati- 
cally to the time evolution of the density matrix. This is 
reasonable since the double occupancy involves processes 
which are generally high in energy and hence fast. Wc 
anticipate that for the cases explored, this approximation 
works well up to energies as large as the Hubbard band 
in the optical conductivity (Fig. |2|) . As the density varies 
the double occupancy shifts from D^ to satisfy the an- 
tiadiabaticity constraint. We define SDit) — D(t) — D^ 
and Spit) and SDit) are linear in /. 

The formal complication of the present approach as 
compared to the standard RPA has its origin in the 
proper adjustment of D to the time evolution of pit), 
i.e. the determination of SDit). This step is achieved 
by expanding the energy functional Eq. (||) up to second 
order in Sp and SD around the saddle point: 



E[p, D]=E + hUp + ^S^LoSp 
+ SDSoSp+ -SD t KoSD 



(4) 



where the bar indicates that we are treating a matrix 
as a column vector and the not indicates evaluation in 
the stationary state. Here we have defined an effective 
one-particle Gutzwiller Hamiltonian fll5|JlC 



- dE 

Opija 



and the matrices 



d 2 E 



dp* jcr dp k i a < 

d 2 E 
dDkdpua 

d 2 E 
dD k dDi ' 



Using the condition of antiadiabaticity 

dE 



dSD 



(5) 

(6) 
(7) 
(8) 

(9) 



in Eq. (^) we obtain a linear relation between Sp and 
SD. Eliminating SD from Eq. (Q) finally yields an ex- 
pansion of the energy as a functional of Sp alone E[p] = 
E\p,D{p)\, 



E[p] =E Q + h 1 Sp + -Sp^Lo ~ SlK^SoW (10) 



This can be regarded as the expansion of an effective 
interacting energy functional in which the interaction 
potential between particles is density dependent. This 
kind of functional often appears in the context of nu- 
clear physics and a well developped machinery exist to 
compute the RPA fluctuations induced by the interac- 
tion. We will only briefly outline here the corresponding 
formalism (for details see Ref. The advantage 

of this method with respect to other methods (eg. dia- 
grammatic) is that the present derivation is solely based 
on the knowledge of an energy functional of a SD density 
matrix which is precisely what the GA provides. 

At the saddle-point h and p can be diagonalized si- 
multaneously. As a result one obtains (/io)fcz = Sutk and 
the density matrix has eigenvalue 1 below the Fermi level 
and eigenvalue above the Fermi level. We will notate 
states below the Fermi level as hole (/i) states and the 
states above the Fermi level as particle states (p). 

Up to linear order the density matrix obeys the equa- 
tion of motion llSjlQl: 



inSp = [h ,5p} + [^.8p,P (0) } 
op 



lf GA ,p {0) ] 



(11) 
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where h is defined as in Eq. (||) but with E instead of E 
(Note that ho = ho). ^-Sp is a short hand notation for 



E 

ph 



dh 



dp. 



ph 



Spph 



c)h 



dp, 



hp 



Sphp 



(12) 



f GA is the GA version of /, i.e. it includes the zo factors 
for intersite matrix elements. 

One can show that particle-particle and hole-hole ma- 
trix elements of Eq. ( p"l| ) are zero and the particle-hole 
(ph) matrix elements of Sp satisfy the well known RPA 
eigenvalue equation. The RPA dynamical matrix can 
be obtained from Eqs. ([To|) , (|Tl|) . Upon diagonalizing 
the RPA matrix by a Bogoliubov transformation one ob- 
tains the eigenvectors — (X^ ) and eigenval- 
ues E^ where the latter correspond to the excitation 
energies of the system. An explicit expression for the 
response functions and a discussion of sum rules can be 
found in Ref. [l5|Jl^| which apply straightforwardly to our 
case. 

The present formalism is well suited for the calcu- 
lation of charge excitations in inhomogcncous doped 
which will be presented elsewhere. 
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systems 

In the following we restrict ourselves to the half-filled 
Hubbard model in the antiferromagnetic Neel state |Q. 
The double occupancy at the RPA level is given by: 
Drpa = J dw£ A (0|n iT |A)<A|n; i |0)5(u;-.E( A )) where the 
integrand is the Lehmann representation of an appropri- 
ately defined density-density correlation function. The 
matrix elements (0|nf|A) for A > can be computed in 
terms of the eigenvectors [ p^[p^ ]. 

In the inset of Fig. [l], we show the GA+RPA double oc- 
cupancy compared with exact results and other approxi- 
mations in a Id system. For small U long range magnetic 
order is not enough to reduce substantially the HF dou- 
ble occupancy from the noninteracting value. RPA on 
top of HF corrects for this but because the starting point 
is quite far from the exact result the correction is not so 
accurate and one gets that HF+RPA overstrikes the ex- 
act double occupancy. On the contrary for the GA only a 
small correction is needed and RPA performs remarkably 
well. Note that U(Drpa — Dhf) is a measure for the 
residual interaction in HF+RPA. In the GA+RPA ap- 
proach such a simple relation is lost but clearly a smaller 
correction of the MF double occupancy suggests a smaller 
residual interaction. 

From the interaction energy UDppA we compute the 
correction to the ground state energy using the coupling 
constant integration trick |2l| . We find very good agree- 
ment with the exact results as shown in Fig. [I]. This holds 
in Id up to intermediate values of U/t and in a 4x4 2d 
cluster for all U jt. The improvement with dimensionality 
is expected as in any MF + RPA computation. 
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FIG. 1. Comparison of the exact ground state energy with 
the various approximate methods discussed in the text for the 
half-filled Hubbard model in Id (upper panel) and 2d (lower 
panel). Exact results in the upper panel are for an infinite 
Id system after Ref. [j^)) whereas approximate results are for 
a 32-site lattice (Finite size errors are estimated to be of the 
order of line width). Exact |^(| and approximated results in 
the lower panel are for a 4x4 cluster. The insets show the 
corresponding curves for the Id double occupancy (exact and 
GA+RPA are almost undistinguishable). 

In order to examine the quality of dynamical corre- 
lation functions we have studied the optical conductiv- 
ity in the GA+RPA approach. As in the HF+RPA 
method |]l5| , ^6|Jl8| the f-sum rule is exactly satisfied with 
the following prescription. The optical conductivity on 
one side of the equality should be computed at the 
GA+RPA (HF+RPA) level and the expectation value 
on the other side (essentially the kinetic energy in our 
case 0) computed at GA (HF) level. 

In this regard the f-sum rule provides also an encourag- 
ing argument that the GA+RPA dynamical correlation 
functions are much more accurate than those obtained 
via the corresponding HF+RPA method. We have com- 
pared the exact kinetic energy of a 4x4 lattice J2f| with 
unrestricted GA and HF results for various hole concen- 
trations and have found that over a wide range of doping 
and on-site correlation U there is almost perfect agree- 
ment between the GA method and exact results. On the 
other hand the HF kinetic energy has an error that for 
example for U = 4t is at least 40 times larger. This is 
not surprising since GA takes into account the correla- 
tion induced reduction of kinetic energy in a much better 



3 



way than HF. 

Fig. H displays a(ui) for a 32-site Hubbard ring and 
half- filling in case of U/t = 4. The onset of excitations 
across the Mott-Hubbard gap is signalled by the appear- 
ance of a large peak in <r(w). We find excellent agreement 
between Monte-Carlo (MC) and GA+RPA whereas the 
excitation energy in HF is clearly overestimated. Note 
that the MC data display an additional hump at approx- 
imately twice the energy of the first peak. Both particle- 
particle scattering processes (not included at RPA level) 
and the failure of the antiadiabaticity assumption for the 
double occupancy at high energies can explain the dis- 
crepancy. We see however that for energies of the order 
of the Hubbard band the method performs very well. 
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FIG. 2. Optical conductivity of a 32-site Hubbard ring 
for U/t = 4 in HF+RPA, GA+RPA and Monte-Carlo (MC) 
after Ref. §e|. The individual peaks of the GA+RPA and 
HF+RPA data have been broadened by a value of 0.25t. 

In conclusion we have presented a time-dependent GA 
for the calculation of dynamical and static quantities in 
the Hubbard model. The approach is conceptually very 
simple and leads to much better agreement with exact 
results than previous approximations. 

As in any computation of fluctuations we are dealing 
with the residual interaction between particles beyond 
the mean-field level. Roughly speaking since the GA con- 
tains ground state correlations not included in a HF wave 
function the residual interaction is a smaller perturbation 
to the MF state and hence it is natural that RPA works 
much better in this case. 

It is interesting to remark that the computation of the 
ground state energy presented here is reminiscent of the 
evaluation of the one in the uniform electron gas based 
on Hubbard type dielectric functions (2^] . Also there the 
computation goes through a density-density correlation 
function and the coupling constant integration trick. Fur- 
thermore the Hubbard local field correction to the dielec- 
tric function takes into account the correlation hole in the 
uniform electron gas whereas the GW projector method 
takes into account similar correlations in the Hubbard 
model. The connection between these two approaches 
deserves further investigation as it may lead to a unified 



approach to strongly correlated systems. 

We greatfully acknowledge many useful discussions 
with C. Di Castro, C. Castellani, M. Grilli and R. Rai- 
mondi. We acknowledge partial financial support from 
INFM. G.S. acknowledges financial support from the 
Deutsche Forschungsgemeinschaft as well as hospitality 
and support from the Dipartimento di Fisica of Univer- 
sita di Roma "La Sapienza" where part of this work was 
carried out. 



[1 
[2 

[3 

[4 

[5, 

ir 

[§ 

[9 
[10 

[11 

[12 

[13 
[14 
[15 

[16 

[17 
[18 
[19 

[20 

[21 

[22' 

[23 

[24 

[25' 

[26 



M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965). 

H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 56, 3582 

(1987). 

W. Metzner and D. Vollhardt, Phys. Rev. Lett. 59, 121 

(1987) . 

W. Metzner and D. Vollhardt, Phys. Rev. B 37, 7382 

(1988) . 

G. Seibold, E. Sigmund, and V. Hizhnyakov, Phys. Rev. 
B 57, 6937 (1998). 

G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 
1362 (1986). 

T. Li, P. Woffle, and P. J. Hirschfeld, Phys. Rev. B 40, 
6817 (1989). 

J. W. Rasul and T. Li, J. Phys. C (Solid St. Phys.) 21, 
5119 (1988). 

M. Lavagna, Phys. Rev. B 41, 142 (1990). 

E. Arrigoni and G. C. Strinati, Phys. Rev. Lett. 71, 3178 

(1993). 

E. Arrigoni and G. C. Strinati, Phys. Rev. B 52, 2428 
(1995). 

R. Raimondi and C. Castellani, Phys. Rev. B 48, R11453 
(1993). 

R. Raimondi, Phys. Rev. B 51, 10154 (1995). 

R. Raimondi, private communication. 

P. Ring and P. Schuck, The nuclear many-body problem 

(Springer- Verlag, New York, 1980). 

J.-P. Blaizot and G. Ripka, Quantum Theory of Finite 
Systems (The MIT Press, Cambridge, Massachusetts, 
1986). 

D. Vollhardt, Rev. Mod. Phys. 56, 99 (1984). 

J. Lorenzana and L. Yu, Phys. Rev. Lett. 70, 861 (1993). 

E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 
(1968). 

G. Fano, F. Ortolani, and A. Parola, Phys. Rev. B 42, 
6877 (1990). 

K. Yonemitsu, A. Bishop, and J. Lorenzana, Phys. Rev. 
B 47, 12059 (1993). 

G. Seibold, C. Castellani, C. D. Castro, and M. Grilli, 
Phys. Rev. B 58, 13506 (1998). 

G. D. Mahan, Many Particle Physics (Plenum, New 
York, 1990). 

B. S. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 
243 (1990). 

E. Dagotto, A. Moreo, F. Ortolani, and J. Riera, Phys. 
Rev. B 45, 10107 (1992). 

N. M. R. Peres, J. M. P. Carmelo, D. K. Campbell, and 
A. W. Sandvik, Z. Phys. B 103, 217 (1997). 



4 



